#This is to get:
#TableA4

library("haven")
library("lfe")
library("rgdal")
library("margins")
library("broom")
library("conleyreg")
library("dplyr")
library("stringr")
library("fwildclusterboot")
library("reshape2")
library("data.table")
library("xtable")


setwd("/Users/bgpopescu/Dropbox/Legacies_Central_Europe/")
setwd("C:/Users/bogdanp/Dropbox/Legacies_Central_Europe/")

#Step0: Table function
make_table<-function(lm_model, keep_vars, no_digits, data){
  #INPUT
  #-lm model
  #-list of vars to keep
  formula <- lm_model$call$formula
  fomula_noxy<-as.character(Reduce(paste, deparse(lm_model$call$formula)))
  fomula_noxy<-gsub(" ", "", fomula_noxy)
  fomula_noxy<-gsub("point_x\\+point_y\\+", "", fomula_noxy)
  fomula_noxy<-as.formula(fomula_noxy)
  new_df<-conleyreg(fomula_noxy, data, 500, lat = "point_x", lon = "point_y")
  df<-tidy(new_df)
  list_df<-list()
  
  for (i in vars_labels){
    iv_name  = i[[1]]
    iv_label = i[[2]]
    df2<-subset(df, term==iv_name)
    list_df[[paste(c("df", iv_name), collapse = '_')]] <- df2}
  
  df2<-data.frame(bind_rows(list_df))
  df2$se<-NA
  df2$se<-ifelse(df2$p.value<0.001,"***",
                 ifelse(df2$p.value<0.01,"**",
                        ifelse(df2$p.value<0.05,"*",
                               ifelse(df2$p.value<0.1,"+",
                                      ""))))
  df2$se_par<-paste("(", round(df2$std.error, no_digits), ")", df2$se, sep="")
  df2<-subset(df2, select = -c(std.error, statistic, p.value, se))
  df2$estimate<-round(df2$estimate, no_digits)
  
  #Bootsrap part
  
  set.seed(1234)
  list_df2<-list()
  for (i in vars_labels){
    iv_name  = i[[1]]
    iv_label = i[[2]]
    a<-Reduce(paste, deparse(lm_model$call$formula))
    #print(iv_name)}
    
    if (grepl(iv_name, a) & iv_name!=word(a, 1)) {
      #print(iv_name)
      boot <- boottest(lm_model, 
                       B = 99999999, 
                       param = iv_name,
                       clustid = "ctrnum",
                       impose_null = F)
      new<-list(summary(boot)$term, as.character(round(summary(boot)$p.value, no_digits)))
      new_df<-data.frame(new)
      names(new_df)<-c("term", "p_value")
      list_df2[[paste(c("df", iv_name), collapse = '_')]] <- new_df}
    else{
      new_df<-data.frame(iv_name, "")
      names(new_df)<-c("term", "p_value")
      list_df2[[paste(c("df", iv_name), collapse = '_')]] <- new_df}}
  
  df3<-data.frame(bind_rows(list_df2))
  
  
  
  df3$p_val_star<-paste("[", as.numeric(df3$p_value), "]",
                        ifelse(as.numeric(df3$p_value)<0.001,"***",
                               ifelse(as.numeric(df3$p_value<0.01),"**",
                                      ifelse(as.numeric(df3$p_value)<0.05,"*",
                                             ifelse(as.numeric(df3$p_value)<0.1,"+",
                                                    "")))), sep = "")
  
  
  
  df3<-subset(df3, select = c(term, p_val_star))
  df4<-left_join(df2, df3, by = c("term" = "term"))
  df5<-melt(setDT(df4), id.vars = c("term"))
  df6 <- df5[order(df5$term),]
  df6$formal<-""
  for (i in vars_labels){
    iv_name  = i[[1]]
    iv_label = i[[2]]
    df6$formal[df6$variable=="estimate" & df6$term==iv_name]<-iv_label
  }
  # 
  no_obs<-data.frame("other", "no_obs", as.character(nobs(lm_model)), "Observations")
  adjr2<-data.frame("other", "adjr2", as.character(round(summary(lm_model)$adj.r.squared, no_digits)), "Adjusted R$^{2}$")
  names(no_obs)<-names(df6)
  names(adjr2)<-names(df6)
  final_list<-list(df6, no_obs, adjr2)
  df7<-do.call("rbind", final_list)
  df7$term_var<-paste(df7$term, df7$variable, sep = "_")
  df8<-subset(df7, select = c("term_var", "formal", "value"))
  return(df8)
}

#Step1: Reading DTA files
nuts_df<-read_dta("./data/nuts_data.dta")

#Step2: Reading shapefile for NUTS polygons (counties)
nuts_sh <- readOGR(dsn="./data/shape_files.gdb",
                   layer="nuts_NONEU_and_EU", use_iconv = TRUE, encoding = "UTF-8")

#Step3: Merge
nuts<-merge(nuts_sh, nuts_df, by.x = "MERGE_ID", by.y = "merge_id")
nuts2<-subset(nuts, otthist==1)
rm(nuts_df, nuts_sh, nuts)

vars_labels<-list(c("p_ottoman", 'Prop. Ottoman'),
                  c("route_density_2", "Trade Route Density"),
                  c("lntmapplic", "Log Trade Mark Certif. per Capita"))
keep_vars2 = c("p_ottoman", "route_density_2", "lntmapplic")



########
#Model1#
########

#Step_M1_1: Remove NAs
mynewspdf2<-nuts2[which(!is.na(nuts2@data$p_ottoman) &
                          !is.na(nuts2@data$lngdpcap)),]

#Step_M1_2: Run regression
reg1_res<-lm(lngdpcap~p_ottoman,
             data=mynewspdf2@data)
summary(reg1_res)

#Step_M1_3: Make Table
sd1<-make_table(reg1_res,
                keep_vars = keep_vars2, 
                no_digits=3, data = mynewspdf2@data)

########
#Model2#
########

#Step_M2_1: Remove NAs
mynewspdf2<-nuts2[which(!is.na(nuts2@data$p_ottoman) &
                          !is.na(nuts2@data$lngdpcap)),]

#Step_M1_2: Run regression
reg2_res<-lm(lngdpcap~p_ottoman + factor(ctrnum),
             data=mynewspdf2@data)
summary(reg2_res)

#Step_M1_3: Make Table
sd2<-make_table(reg2_res,
                keep_vars = keep_vars2, 
                no_digits=3, data = mynewspdf2@data)


########
#Model3#
########

#Step_M3_1: Remove NAs
mynewspdf2<-nuts2[which(!is.na(nuts2@data$lngdpcap) &
                          !is.na(nuts2@data$p_ottoman) &
                          !is.na(nuts2@data$ctrnum) &
                          !is.na(nuts2@data$tavg) &
                          !is.na(nuts2@data$lograin) &
                          !is.na(nuts2@data$logelev) &
                          !is.na(nuts2@data$logslope) & 
                          !is.na(nuts2@data$wheat_rainfed) &
                          !is.na(nuts2@data$log_sea_lines_distance) & 
                          !is.na(nuts2@data$log_large_rivers_distance) & 
                          !is.na(nuts2@data$point_x) & 
                          !is.na(nuts2@data$point_y)),]


#Step_M1_2: Run regression
reg3_res<-lm(lngdpcap~p_ottoman +
               tavg + lograin + logelev +
               logslope + wheat_rainfed + log_sea_lines_distance +
               log_large_rivers_distance + point_x + point_y +
               factor(ctrnum),
             data=mynewspdf2@data)
summary(reg3_res)

#Step_M1_3: Make Table
sd3<-make_table(reg3_res,
                keep_vars = keep_vars2, 
                no_digits=3, data = mynewspdf2@data)
sd3

#########
#Model 4#
#########

#Step_M4_1: Remove NAs
mynewspdf2<-nuts2[which(!is.na(nuts2@data$lngdpcap) &
                          !is.na(nuts2@data$p_ottoman) &
                          !is.na(nuts2@data$route_density_2) &
                          !is.na(nuts2@data$ctrnum) &
                          !is.na(nuts2@data$tavg) &
                          !is.na(nuts2@data$lograin) &
                          !is.na(nuts2@data$logelev) &
                          !is.na(nuts2@data$logslope) & 
                          !is.na(nuts2@data$wheat_rainfed) &
                          !is.na(nuts2@data$log_sea_lines_distance) & 
                          !is.na(nuts2@data$log_large_rivers_distance) & 
                          !is.na(nuts2@data$point_x) & 
                          !is.na(nuts2@data$point_y)),]

#Step_M1_2: Run regression
#Note that I take out lat lon. Conley Reg does not work with lat lon inside
reg4_res<-lm(lngdpcap~p_ottoman + route_density_2 +
               tavg + lograin + logelev +
               logslope + wheat_rainfed + log_sea_lines_distance +
               log_large_rivers_distance + 
               point_x + point_y + 
               factor(ctrnum),
             data=mynewspdf2@data)
summary(reg4_res)

#Step_M1_3: Make Table
sd4<-make_table(reg4_res,
                keep_vars = keep_vars2, 
                no_digits=3, data = mynewspdf2@data)
sd4


#########
#Model 8#
#########


#Step_M8_1: Remove NAs
mynewspdf2<-nuts2[which(!is.na(nuts2@data$firmscap) &
                          !is.na(nuts2@data$p_ottoman)  &
                          !is.na(nuts2@data$ctrnum) &
                          !is.na(nuts2@data$tavg) &
                          !is.na(nuts2@data$lograin) &
                          !is.na(nuts2@data$logelev) &
                          !is.na(nuts2@data$logslope) & 
                          !is.na(nuts2@data$wheat_rainfed) &
                          !is.na(nuts2@data$log_sea_lines_distance) & 
                          !is.na(nuts2@data$log_large_rivers_distance) & 
                          !is.na(nuts2@data$point_x) & 
                          !is.na(nuts2@data$point_y)),]


#Step_M8_2: Run regression
reg8_res<-lm(firmscap~p_ottoman + 
               tavg + lograin + logelev +
               logslope + wheat_rainfed + log_sea_lines_distance +
               log_large_rivers_distance + 
               point_x + point_y + 
               factor(ctrnum),
             data=mynewspdf2@data)
summary(reg8_res)

#Step_M8_3: Make Table
sd8<-make_table(reg8_res,
                keep_vars = keep_vars2, 
                no_digits=3, data = mynewspdf2@data)
sd8

#########
#Model 9#
#########

#Step_M9_1: Remove NAs
mynewspdf2<-nuts2[which(!is.na(nuts2@data$firms10cap) &
                          !is.na(nuts2@data$p_ottoman)  &
                          !is.na(nuts2@data$ctrnum) &
                          !is.na(nuts2@data$tavg) &
                          !is.na(nuts2@data$lograin) &
                          !is.na(nuts2@data$logelev) &
                          !is.na(nuts2@data$logslope) & 
                          !is.na(nuts2@data$wheat_rainfed) &
                          !is.na(nuts2@data$log_sea_lines_distance) & 
                          !is.na(nuts2@data$log_large_rivers_distance) & 
                          !is.na(nuts2@data$point_x) & 
                          !is.na(nuts2@data$point_y)),]

#Step_M9_2: Run regression
reg9_res<-lm(firms10cap~p_ottoman + 
               tavg + lograin + logelev +
               logslope + wheat_rainfed + log_sea_lines_distance +
               log_large_rivers_distance + 
               point_x + point_y + 
               factor(ctrnum),
             data=mynewspdf2@data)
summary(reg9_res)

#Step_M9_3: Make Table
sd9<-make_table(reg9_res,
                keep_vars = keep_vars2, 
                no_digits=3, data = mynewspdf2@data)
sd9


###########
#Model 10#
##########

#Step_M10_1: Remove NAs
mynewspdf2<-nuts2[which(!is.na(nuts2@data$lntmapplic) &
                          !is.na(nuts2@data$p_ottoman)  &
                          !is.na(nuts2@data$ctrnum) &
                          !is.na(nuts2@data$tavg) &
                          !is.na(nuts2@data$lograin) &
                          !is.na(nuts2@data$logelev) &
                          !is.na(nuts2@data$logslope) & 
                          !is.na(nuts2@data$wheat_rainfed) &
                          !is.na(nuts2@data$log_sea_lines_distance) & 
                          !is.na(nuts2@data$log_large_rivers_distance) & 
                          !is.na(nuts2@data$point_x) & 
                          !is.na(nuts2@data$point_y)),]

#Step_M10_2: Run regression
reg10_res<-lm(lntmapplic~p_ottoman + 
                tavg + lograin + logelev +
                logslope + wheat_rainfed + log_sea_lines_distance +
                log_large_rivers_distance + 
                point_x + point_y + 
                factor(ctrnum),
              data=mynewspdf2@data)
summary(reg10_res)

#Step_M10_3: Make Table
sd10<-make_table(reg10_res,
                keep_vars = keep_vars2, 
                no_digits=3, data = mynewspdf2@data)
sd10



##########
#Model 11#
##########

#Step_M11_1: Remove NAs
mynewspdf2<-nuts2[which(!is.na(nuts2@data$lngdpcap) &
                          !is.na(nuts2@data$lntmapplic) &
                          !is.na(nuts2@data$p_ottoman)  &
                          !is.na(nuts2@data$ctrnum) &
                          !is.na(nuts2@data$tavg) &
                          !is.na(nuts2@data$lograin) &
                          !is.na(nuts2@data$logelev) &
                          !is.na(nuts2@data$logslope) & 
                          !is.na(nuts2@data$wheat_rainfed) &
                          !is.na(nuts2@data$log_sea_lines_distance) & 
                          !is.na(nuts2@data$log_large_rivers_distance) & 
                          !is.na(nuts2@data$point_x) & 
                          !is.na(nuts2@data$point_y)),]

#Step_M11_2: Run regression
reg11_res<-lm(lngdpcap~p_ottoman + lntmapplic+
                tavg + lograin + logelev +
                logslope + wheat_rainfed + log_sea_lines_distance +
                log_large_rivers_distance + 
                point_x + point_y + 
                factor(ctrnum),
              data=mynewspdf2@data)
summary(reg11_res)

#Step_M11_3: Make Table
sd11<-make_table(reg11_res,
                 keep_vars = keep_vars2, 
                 no_digits=3, data = mynewspdf2@data)
sd11

################
#MODEL 12 FILES#
################

#Step1: Reading DTA files
educ_df<-read_dta("./data/educ_data.dta")

#Step2: Reading shapefile for NUTS polygons (counties)
nuts_sh <- readOGR(dsn="./data/shape_files.gdb",
                   layer="nuts_EU_2001", use_iconv = TRUE, encoding = "UTF-8")
#Step3: Merge
nuts_edu<-merge(nuts_sh, educ_df, by.x = "MERGE_ID_EDU", by.y = "MERGE_ID_EDU")


##########
#MODEL 12#
##########
#Step_M12_1: Remove NAs
mynewspdf_2<-nuts_edu[which(!is.na(nuts_edu@data$secondaryed2) &
                              !is.na(nuts_edu@data$p_ottoman) &
                              !is.na(nuts_edu@data$route_density_2) &
                              !is.na(nuts_edu@data$ctrnum) &
                              !is.na(nuts_edu@data$tavg) &
                              !is.na(nuts_edu@data$lograin) &
                              !is.na(nuts_edu@data$logelev) &
                              !is.na(nuts_edu@data$logslope) & 
                              !is.na(nuts_edu@data$point_x) & 
                              !is.na(nuts_edu@data$point_y) &
                              !is.na(nuts_edu@data$wheat_rainfed) &
                              !is.na(nuts_edu@data$log_sea_lines_distance) &
                              !is.na(nuts_edu@data$log_large_rivers_distance)),]

#Step_M12_2: Run regression
reg12_res<-lm(secondaryed2~p_ottoman +
                tavg + lograin + logelev +
                logslope + 
                wheat_rainfed + log_sea_lines_distance +
                log_large_rivers_distance + 
                point_x + point_y + 
                factor(ctrnum),
              data=mynewspdf_2@data)
summary(reg12_res)

#Step_M12_3: Make Table
sd12<-make_table(reg12_res,
                 keep_vars = keep_vars2, 
                 no_digits=3, data = mynewspdf_2@data)
sd12

######################
#Organizing the Table#
######################

#Step1: Collecting the models
new<-c(sd2$term_var, sd3$term_var, sd4$term_var,
       sd8$term_var, sd9$term_var,sd10$term_var,
       sd11$term_var,sd12$term_var)

#Step2: Ordering rows
new2<-unique(new)
new_df<-data.frame(new2)
new_df$order[new_df$new2=="p_ottoman_estimate"]<-1
new_df$order[new_df$new2=="p_ottoman_se_par"]<-2
new_df$order[new_df$new2=="p_ottoman_p_val_star"]<-3
new_df$order[new_df$new2=="route_density_2_estimate"]<-4
new_df$order[new_df$new2=="route_density_2_se_par"]<-5
new_df$order[new_df$new2=="route_density_2_p_val_star"]<-6
new_df$order[new_df$new2=="lntmapplic_estimate"]<-7
new_df$order[new_df$new2=="lntmapplic_se_par"]<-8
new_df$order[new_df$new2=="lntmapplic_p_val_star"]<-9
new_df$order[new_df$new2=="other_no_obs"]<-17
new_df$order[new_df$new2=="other_adjr2"]<-18
new_df <- new_df[order(new_df$order),]
names(new_df)[names(new_df) == "new2"] <- "term_var"

#Step3: Subsetting dataframes
sd1<-subset(sd1, select = c("term_var", "value"))
names(sd1)[names(sd1) == "value"] <- "model_01"
sd2<-subset(sd2, select = c("term_var", "value"))
names(sd2)[names(sd2) == "value"] <- "model_02"
sd3<-subset(sd3, select = c("term_var", "value"))
names(sd3)[names(sd3) == "value"] <- "model_03"
sd4<-subset(sd4, select = c("term_var", "value"))
names(sd4)[names(sd4) == "value"] <- "model_04"
sd8<-subset(sd8, select = c("term_var", "value"))
names(sd8)[names(sd8) == "value"] <- "model_08"
sd9<-subset(sd9, select = c("term_var", "value"))
names(sd9)[names(sd9) == "value"] <- "model_09"
sd10<-subset(sd10, select = c("term_var", "value"))
names(sd10)[names(sd10) == "value"] <- "model_10"
sd11<-subset(sd11, select = c("term_var", "value"))
names(sd11)[names(sd11) == "value"] <- "model_11"
sd12<-subset(sd12, select = c("term_var", "value"))
names(sd12)[names(sd12) == "value"] <- "model_12"

#Step4 Joining the dataframes
final<-left_join(new_df, sd1, by='term_var') %>%
  left_join(., sd2, by='term_var') %>%
  left_join(., sd3, by='term_var') %>%
  left_join(., sd4, by='term_var') %>%
  left_join(., sd8, by='term_var') %>%
  left_join(., sd9, by='term_var') %>%
  left_join(., sd10, by='term_var') %>%
  left_join(., sd11, by='term_var') %>%
  left_join(., sd12, by='term_var')

#Step5 Giving formal names
final$formal<-NA
final$formal[final$term_var=="p_ottoman_estimate"]<-"Prop. Ottoman"
final$formal[final$term_var=="route_density_2_estimate"]<-"Trade Route Density"
final$formal[final$term_var=="lntmapplic_estimate"]<-"Log Trade Mark Certif. per Capita"
final$formal[final$term_var=="other_no_obs"]<-"Observations"
final$formal[final$term_var=="other_adjr2"]<-"Ajusted $R^2$"
new<-c("cfe", 12, "No", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Country FE")
final2<-rbind(final, new)
new2<-c("model", 13, rep("OLS", 9), "Model")
final3<-rbind(final2, new2)
new3<-c("ctrls", 14, "No", "No", rep("Geo+Pos", 7), "Controls")
final4<-rbind(final3, new3)
new4<-c("years", 15, rep("All", 9), "Years")
final5<-rbind(final4, new4)

#Step6 Merging the dataframes
final_tab<-subset(final5, select = c(formal, model_01, model_02, model_03, 
                                     model_04, model_08, model_09, model_10, 
                                     model_11, model_12))

#########################
#MARGINAL EFFECTS MODELS#
#########################

#Step_M5_1: Run regression
reg5_res<-lm(lngdpcap~p_ottoman * yearpos+ 
               tavg + lograin + logelev +
               logslope + wheat_rainfed + log_sea_lines_distance +
               log_large_rivers_distance + point_x + point_y +
               factor(ctrnum),
             data=mynewspdf2@data)
summary(reg5_res)

#Step_M5_2: Saving coefficients and SE
coef_reg5_res<-summary(margins(reg5_res,  vce = "delta", variables = "p_ottoman", at = list(yearpos = 1650)))$AME
se_reg5_res<-summary(margins(reg5_res,  vce = "delta", variables = "p_ottoman", at = list(yearpos = 1650)))$SE

#Step M5_3: Bootstrap SE - 
#No implementation for bootstrap margins in R
#Run in Stata the following code
#egen ctrnum1 = group(ctrnum)
#reg lngdpcap c.p_ottoman##c.yearpos tavg lograin logelev logslope wheat_rainfed log_sea_lines_distance log_large_rivers_distance point_x point_y i.ctrnum1, cluster(ctrnum1)
#margins, dydx(p_ottoman) at(yearpos = 1650)
#boottest, margins  cluster(ctrnum1) nonull nograph seed(1234)
reg5_res_f_pvalue<-0.0180
reg5_res_pvalue_f_star<-paste("[", round(reg5_res_f_pvalue, digits = 3),
                              ifelse(reg5_res_f_pvalue<0.001,"]***",
                                     ifelse(reg5_res_f_pvalue<0.01,"]**",
                                            ifelse(reg5_res_f_pvalue<0.05,"]*",
                                                   ifelse(reg5_res_f_pvalue<0.1,"]+",
                                                          "]")))), sep="")

#Step M6_1: Run regression
reg6_res<-lm(lngdpcap~p_ottoman * yearpos+ 
               tavg + lograin + logelev +
               logslope + wheat_rainfed + log_sea_lines_distance +
               log_large_rivers_distance + point_x + point_y +
               factor(ctrnum),
             data=mynewspdf2@data)
summary(reg6_res)

#Step_M6_2: Saving coefficients and SE
coef_reg6_res<-summary(margins(reg6_res,  vce = "delta", variables = "p_ottoman", at = list(yearpos = 1700)))$AME
se_reg6_res<-summary(margins(reg6_res,  vce = "delta", variables = "p_ottoman", at = list(yearpos = 1700)))$SE


#Step M6_3: Bootstrap SE
#No implementation for bootstrap margins in R
#Run in Stata the following code
#reg lngdpcap c.p_ottoman##c.yearpos tavg lograin logelev logslope wheat_rainfed log_sea_lines_distance log_large_rivers_distance point_x point_y i.ctrnum if (ctr != "TR" & otthist==1), cluster(ctrnum)
#margins, dydx(p_ottoman) at(yearpos = 1700)
#boottest, margins  cluster(ctrnum) nonull nograph reps (9999) seed(1234)
reg6_res_f_pvalue<-0.0000

reg6_res_pvalue_f_star<-paste("[", round(reg6_res_f_pvalue, digits = 3),
                              ifelse(reg6_res_f_pvalue<0.001,"]***",
                                     ifelse(reg6_res_f_pvalue<0.01,"]**",
                                            ifelse(reg6_res_f_pvalue<0.05,"]*",
                                                   ifelse(reg6_res_f_pvalue<0.1,"]+",
                                                   "]")))), sep="")


#Step M7_1: Run regression
reg7_res<-lm(lngdpcap~p_ottoman * yearpos+ 
               tavg + lograin + logelev +
               logslope + wheat_rainfed + log_sea_lines_distance +
               log_large_rivers_distance + point_x + point_y +
               factor(ctrnum),
             data=mynewspdf2@data)
summary(reg7_res)

#Step M7_2: Saving SE
coef_reg7_res<-summary(margins(reg7_res,  vce = "delta", variables = "p_ottoman", at = list(yearpos = 1750)))$AME
se_reg7_res<-summary(margins(reg7_res,  vce = "delta", variables = "p_ottoman", at = list(yearpos = 1750)))$SE

#Step M7_3: Bootstrap SE
#No implementation for bootstrap margins in R
#Run in Stata the following code
#reg lngdpcap c.p_ottoman##c.yearpos tavg lograin logelev logslope wheat_rainfed log_sea_lines_distance log_large_rivers_distance point_x point_y i.ctrnum1 if (otthist==1), cluster(ctrnum1)
#margins, dydx(p_ottoman) at(yearpos = 1750)
#boottest, margins  cluster(ctrnum1) nonull nograph reps (9999) seed(1234)

reg7_res_f_pvalue<-0.3564

reg7_res_pvalue_f_star<-paste("[", round(reg7_res_f_pvalue, digits = 3),
                              ifelse(reg7_res_f_pvalue<0.001,"]***",
                                     ifelse(reg7_res_f_pvalue<0.01,"]**",
                                            ifelse(reg7_res_f_pvalue<0.05,"]*",
                                                   ifelse(reg7_res_f_pvalue<0.1,"]*",
                                                   "]")))), sep="")
######################
#Organizing the Table#
######################
#Step1 Making lists for marginal effects
new_col_1650<-list(round(coef_reg5_res, 3), 
                   "",
                   reg5_res_pvalue_f_star,
                   "", "", "", 
                   "", "", "", 
                   nobs(reg5_res), "", "Yes",
                   "OLS", "Geo + Pos", "1650")

new_col_1700<-list(round(coef_reg6_res, 3),
                   "",
                   reg6_res_pvalue_f_star,
                   "", "", "", 
                   "", "", "", 
                   nobs(reg6_res), "","Yes",
                   "OLS", "Geo + Pos", "1700")

new_col_1750<-list(round(coef_reg7_res, 3), 
                   "",
                   reg7_res_pvalue_f_star,
                   "", "", "", 
                   "", "", "", 
                   nobs(reg7_res), "","Yes", 
                   "OLS", "Geo + Pos", "1750")

#Step2 Adding the lists to the dataframe
final_tab$new_col_1650<-new_col_1650
names(final_tab)[names(final_tab) == "new_col_1650"]<-"model_05"
final_tab$new_col_1700<-new_col_1700
names(final_tab)[names(final_tab) == "new_col_1700"]<-"model_06"
final_tab$new_col_1750<-new_col_1750
names(final_tab)[names(final_tab) == "new_col_1750"]<-"model_07"

#Step3 Order columns
final1<-final_tab[ , order(names(final_tab))]
new_row<-names(final1)

#Step4 Create numbers for columns (stargazer style)
clean_new_row<-c()
for (i in new_row){
  x2<-as.numeric(gsub("model_", "", i))
  x3<-as.character(x2)
  print(x3)
  if (is.na(x3)){
    x3<-""
    clean_new_row<-append(clean_new_row, x3)}
  else{
    x4<-paste("(", x3, ")", sep="")
    clean_new_row<-append(clean_new_row, x4)}}
clean_new_row[clean_new_row!= "(NA)"]

#Step5 Binding row numbers to dataframe
final2<-rbind(clean_new_row, final1)

#Step6 Giving DV names
names(final2)<-c("", "\\shortstack{Log\\\\GDP/cap}",
              rep("\\shortstack{Log\\\\GDP/cap}",6),
              "\\shortstack{Firms\\\\per\\\\Capita}",
              "\\shortstack{Firms\\\\> 10 employees\\\\per Capita}",
              "\\shortstack{Trade Mark\\\\Certifications\\\\per Capita}",
              "\\shortstack{Log\\\\GDP/cap}",
              "\\shortstack{Share. Pop.\\\\with Sec.\\\\Education}")


#Step7 centering coefficients
object_latex<-xtable(final2, type = "latex", 
                     caption = "Results",
                     #Note there are X columns to be aligned below.
                     #This is because of the index column.
                     align = "llcccccccccccc")

comment <- list(pos = list(0), command = NULL)
comment$pos[[1]] <- c(nrow(object_latex))
comment$command <- c(paste("\\bottomrule\n",
                           "\\multicolumn{13}{l} {\\parbox[t]{37cm}\n",
                           "{\\footnotesize \\textit{Note}: 
                            All columns present the marginal effects of the proportion Ottoman indicator (range 0-1) 
                            on ln(GDP/capita) in 2016 or on the other dependent variables at the regional level. 
                            Conley standard errors using a radius of allowed spatial dependence of 500km
                            are reported in parentheses.
                            Wild Cluster Bootstrapping p-values are reported in square brackets.
                            Geo controls include average temperature, ln(average altitude), ln(average slope), 
                            ln(average precipitation), wheat growth suitability, distance from rivers, 
                            and distance from the sea. Pos controls include latitude and longitude. 
                            Models 5, 6, and 7 present results from models in which proportion Ottoman is 
                            interacted with the average year of Ottoman rule, and marginal effects in 
                            1650, 1700, and 1750 are calculated. *** = p<.001; ** = p<.01, *=p<.05, +=p<.10.
                           *** = p<.001; ** = p<.01, *=p<.05, +=p<.10.}} \\\\ \n", sep = ""))

#Step8 Printing latex tex
print(object_latex,
      #The following line controls where the lines are drawn
      hline.after=c(-1, 1, 10),
      floating = FALSE,
      file = "./Paper/tables/tableA4.tex", 
      #caption.placement = "top",
      add.to.row = comment,
      booktabs = TRUE,
      include.rownames=FALSE,
      sanitize.text.function=function(x){x})